General remarks

The R pacakge enerscape computes the energy landscape based on models predicting the energy costs of transport. Two models are currently implemented, the ARC model for terrestrial, legged animals from (Pontzer 2016) and the model for cycling expenditures from (Prampero et al. 1979). Adding other models is relatively easy and can be done by writing a custom function in the enerscape_internals.R script. Currently, there is no implementation for user-specified models, but I plan to add this in the future.

An extensive derivation of the ARC model can be found in Pontzer (2016). Here, it is sufficient to say that the ARC model accounts for the energy spent for cross-bridge cycling, generating tension, as well as associated activation-relaxation processes that restore cellular ion gradients - mostly Ca\(^{2+}\) necessary for cross-bridge cycling. Notably, predictions of the ARC model fit very well observed data. In enerscape, the dataset from Pontzer (2016) can be access as the variable pontzer:

library(tidyverse)
library(raster)
library(sf)
library(enerscape)

#' Function to plot the areas
#' 
#' @param x is a raster
#' @param poly is a shapefile polygon
plot_area <- function(x, poly, mask = NULL, col = NULL, void = FALSE) {
  if (is.null(col)) {
    plot(x, axes = !void, box = !void)
  } else {
    plot(x, col = col, axes = !void, box = !void)
  }
  if (!is.null(mask)) {
    plot(mask, add = T, legend = FALSE,
         col = adjustcolor("white", alpha.f = 0.5), 
         axes = !void, 
         box = !void)
  }
  plot(poly$geometry, add = TRUE)
}

#' ARC model predictions
#'
#' @param mass body mass of the animal (kg)
#' @param slope incline of the terrain (degree)
ARC <- function(mass, slope) {
  E_ar <- 8 * mass ^ (-0.34)
  E_mec <- 100 * (1 + sin((2 * slope - 74) / 180 * pi)) * mass ^ (-0.12)
  E <- (E_ar + E_mec)
  return(E)
}

pontzer %>% 
  as_tibble %>% 
  mutate(ARC = ARC(Mass, Incline)) %>% 
  ggplot() +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed") +
  geom_point(aes(Cost.of.Transport, ARC, col = Incline)) +
  scale_color_gradient(low = "steelblue", high = "tomato") +
  scale_x_log10() +
  scale_y_log10() +
  xlab(expression("E"["COT"]*" (J kg"^"-1"*"m"^"-1"*")")) +
  ylab("ARC predicitons") +
  theme_bw()

Computing energy landscapes

To compute energy landscapes using enerscape, we need an elevation rasterLayer, normally called a digital elevation model (DEM). The raster should have equal resolution on both horizontal and vertical axes and a planar coordinate reference systems in meters. The following is a raster that have these characteristics:

dem <- raster("../data/dem-vignette.tif")
res(dem)

## [1] 100 100

crs(dem)

## CRS arguments:
##  +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs

We can then compute the energy landscape of the area for a general animal of 50 kg using the enerscape() function:

en <- enerscape(dem, 50)

##  - Raster cells are assumed to have same horizontal and vertical resolution and with planar coordinate reference system (e.g. UTM)

##   | Calculating slope

##   | Calculating work

##   | Calculating conductance (1 / work)

##  - Do not use slope with negative values for distance calculations

The output is an enerscape object, which is a list with elements:

• The neighbors used to compute the transition matrix among cells, either 4 (in chess, rook’s move of 1), 8 (king’s move), or 16 (king’s + knigth’s move).
• The body mass of the animal (kg).
• A rasterStack with DEM, Slope, Work (cost of travel), and Conductance (1 / Work).
• A conductance transition matrix that summarize the conductance of moving between neighboring cells.

The function en_extrapolation() can be used to highlight where ARC predictions extrapolate from the testing data. A rasterLayer masking extrapolations for the incline is returned if available:

extr <- en_extrapolation(en, plot = FALSE)
extr

## $`Mass extrapolated`
## [1] FALSE
## 
## $`Slope extrapolated`
## [1] TRUE
## 
## $`Slope extrapolation`
## class      : RasterLayer 
## dimensions : 367, 509, 186803  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 852140, 903040, 4659980, 4696680  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : Slope 
## values     : 1, 1  (min, max)

plot(dem)
plot(extr$`Slope extrapolation`, 
     add = TRUE, 
     col = adjustcolor("blue", alpha.f = 0.3),
     legend = FALSE)
arrows(888523, 4674275, 878791, 4671549)

The algorithm to compute least-cost paths from the R package gdistance is implemented in the en_lcp() functio0n (Etten 2017). en_lcp() takes as input two points for which the least-cost path need to be evaluated. Alternatively, random points are generated if simulate_random_points = TRUE:

lcp <- enerscape::en_lcp(en, simulate_random_points = TRUE)

## 1 of 10 random points ...

## 2 of 10 random points ...

## 3 of 10 random points ...

## 4 of 10 random points ...

## 5 of 10 random points ...

## 6 of 10 random points ...

## 7 of 10 random points ...

## 8 of 10 random points ...

## 9 of 10 random points ...

## 10 of 10 random points ...

The random walk algorithm from gdistance is also implemented as en_passage():

walk <- enerscape::en_passage(en,
                              simulate_random_points = TRUE,
                              theta = 4)

## 1 of 10random points ...

## 2 of 10random points ...

## 3 of 10random points ...

## 4 of 10random points ...

## 5 of 10random points ...

## 6 of 10random points ...

## 7 of 10random points ...

## 8 of 10random points ...

## 9 of 10random points ...

## 10 of 10random points ...

walk

## $Origins
##            x       y
##  [1,] 861890 4687030
##  [2,] 868490 4690430
##  [3,] 852890 4682130
##  [4,] 894290 4665930
##  [5,] 855990 4678530
##  [6,] 866890 4686230
##  [7,] 900490 4692930
##  [8,] 893790 4681730
##  [9,] 874190 4666630
## [10,] 857690 4681530
## 
## $Destinations
##            x       y
##  [1,] 902690 4690630
##  [2,] 852490 4695830
##  [3,] 863090 4677930
##  [4,] 864190 4674430
##  [5,] 895690 4679530
##  [6,] 869790 4667330
##  [7,] 889990 4686530
##  [8,] 902590 4666030
##  [9,] 855890 4679330
## [10,] 863790 4662530
## 
## $Passages
## class      : RasterStack 
## dimensions : 367, 509, 186803, 10  (nrow, ncol, ncell, nlayers)
## resolution : 100, 100  (x, y)
## extent     : 852140, 903040, 4659980, 4696680  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## names      : layer.1, layer.2, layer.3, layer.4, layer.5, layer.6, layer.7, layer.8, layer.9, layer.10 
## min values :       0,       0,       0,       0,       0,       0,       0,       0,       0,        0 
## max values :       0,       0,       0,       0,       0,       0,       0,       0,       0,        0 
## 
## 
## $`Cumulative net passage`
## class      : RasterLayer 
## dimensions : 367, 509, 186803  (nrow, ncol, ncell)
## resolution : 100, 100  (x, y)
## extent     : 852140, 903040, 4659980, 4696680  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=32 +datum=WGS84 +units=m +no_defs 
## source     : memory
## names      : layer 
## values     : NA, NA  (min, max)

As this may be quite slow, I implemented another way of computing the overall connectivity between two points - or for the whole landscape - using Circuitscape and Omniscape. The initialization files for these to algorithms can be written using circuitscape_skeleton() and omniscape_skeleton(), which create circuitscape.ini and omniscape.ini to be launched from within a Julia terminal. See the example [Landscape connectivity for the Marsican bear (Ursus arctos marsicanus) in the Sirente-Velino Regional Park: Circuitscape and Omniscape implementations.] for details.